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Abstract The correct estimation of thermal properties 
of ultra-scaled CMOS and thermoelectric semiconduc- 
tor devices demands for accurate phonon modeling in 
such structures. This work provides a detailed descrip- 
tion of the modified valence force field (MVFF) method 
to obtain the phonon dispersion in zinc-blende semi- 
conductors. The model is extended from bulk to nano- 
wires after incorporating proper boundary conditions. 
The computational demands by the phonon calculation 
increase rapidly as the wire cross-section size increases. 
It is shown that the nanowire phonon spectrum differ 
considerably from the bulk dispersions. This manifests 
itself in the form of different physical and thermal prop- 
erties in these wires. We believe that this model and 
approach will prove beneficial in the understanding of 
the lattice dynamics in the next generation ultra-scaled 
semiconductor devices. 

Keywords Dynamical matrix • Nanowire • Phonons • 
Valence Force Field 

PACS 63.20.-e Phonons in crystal lattices • 63.22.Gh 
Nanotubes and nanowires 



1 Introduction 

The lattice vibration modes known as 'phonons' de- 
termine many important properties in semiconductors 
like, (i) the phonon limited low field carrier mobility in 
mosFETs [HE], (ii) the lattice thermal conductivity in 
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semiconductors which plays an important role in ther- 
moelectric design [3U21IH| , and (iii) the structural stabil- 
ity of ultra-thin semiconductor nanowires [6] . A physics- 
based method to calculate the phonon dispersion in 
semiconductors is required to understand and link all 
these issues. As the device size approach the nanometer 
scale and as the number of atoms in the structure be- 
come countably finite, a continuum material description 
is no longer accurate. This work provides a complete 
and elaborate description of an atomistic phonon cal- 
culation method based on 'Valence Force Field' (VFF) 
model [71IHUS] , a frozen phonon approach, with applica- 
tion to bulk and nanowire structures. 

A variety of methods have been reported in the liter- 
ature for the calculation of the phonon spectrum such as 
the Valence Force Field (VFF) method and its variants 
[MUTED], Bond Charge Model (BCM) [HUH], Density 
Functional Method [S1[T3"]. etc. We focus on VFF meth- 
ods in this work. There are multiple reasons for using a 
VFF based model: (a) in covalent bonded crystals, like 
Si, Ge, GaAs, simple VFF potentials are sufficient to 
match the experimental data jTUj, (b) valence coordi- 
nates and hence the potential energy (U) depend only 
on the relative positions of the atoms and are inde- 
pendent of rigid translations and rotations of the solid, 
and (c) it is easy to extend the model to confined ultra- 
scaled structures made of few atoms since the interac- 
tions are at the atomic level. 

The original Keating VFF model [7] describe the 
LA, LO and TO phonons reasonably well in zinc-blende 
materials, however, it does not produce the flatness in 
the TA branch in Si, Ge [THIS]. Also the limitation of 
the model to correctly describe the elastic constants 
(Cll, C12, C44) in these materials also limits its use 
|S]. In order to extend the available VFF models to 
nanostructures, we need to identify the models which 
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can correctly describe the phonons in the entire Bril- 
louin zone (BZ) in zinc-blende semiconductors. There 
are older works where as many as six parameters |14) 
have been used in VFF to obtain the correct phonon 
dispersions. However, the task of this work has been 
to obtain a VFF model which (i) can capture the cor- 
rect physics and (ii) is computationally not very ex- 
pensive. Hence such a model can be extended to nano- 
structurcs like nanowires, ultra-thin-bodies, etc. To this 
end we have identified two VFF models which satisy the 
requirements. The modified VFF (MVFF) model pre- 
sented here combines these two following models, (i) 
VFF model from Sui et. al [5] which is suitable for non- 
polar materials like Si and Ge and (ii) VFF model from 
Zunger et. al [9] which is suitable for treating polar 
materials like GaP, GaAs, etc. This extended model is 
called the 'MVFF model' in this study. 

The main focus of this work is to show the imple- 
mentation of VFF models for phonon calculation in 
zinc-blende (diamond) lattices. We present the details 
on the atomic groups which make up the interactions, 
the application of boundary conditions in the nano- 
structures, the eigen value problem, the computational 
requirements and the evaluation of lattice properties. 
We benchmarked the model for variety of zinc-blende 
materials like Si, Ge, GaP, GaAs, etc. In this paper we 
present the results using Si (sometimes Ge too) as a 
specific example. We also present a comparison of the 
Keating VFF model [7j with the present MVFF model 
for Si to elucidate the differences in physical results and 
their computational requirements. 

Previous theoretical works have reported the cal- 
culation of phonon dispersions in SiNWs using a con- 
tinuum elastic model and Boltzmann transport equa- 
tion [T5], atomistic first principle methods like DFPT 
(Density Functional Perturbation Theory) |7jIiT3"lil6lll7| 
and atomistic frozen phonon approaches like Keating- 
VFF (KVFF) (HUH]. Thermal conductivity in SiNWs 
has been studied previously using the KVFF model [4] 

Eg. 

This paper has been arranged in the following sec- 
tions. The MVFF theory (a 'frozen phonon' method) 
for the phonon dispersion in zinc-blende semiconduc- 
tors is reviewed in Sec. [2] It provides details about 
the total potential energy (U) of the crystals in the 
MVFF model (Sec. |2.1[ ), construction of the dynam- 
ical matrix (DM) (Sec. 2.2|, application of boundary 



conditions to the DM (Sec. |2.3|), so lution of the result- 
ing eigen value problem (Sec. 2.4 1, and calculation of 
sound velocity (V sn d) (Sec. |2.5[ ), lattice thermal conduc- 
tance (<7/) (Sec. 2.6) and the mode Griineisen parame- 



ters (Sec. |2.7[ ) using the phonon spectrum. The compu- 
tational details for the calculation of phonon dispersion 



are presented in Sec. [3] Different aspects of the dynami- 
cal matrix like the size, fill-factor, sparsity pattern, etc., 
are provided in Sec. |3.1| Timing analysis for the assem- 
bly of the DM are shown in Sec |3.2| Section [4] presents, 
a benchmark of the MVFF results against experimen- 
tal data for different semiconductors (Sec. 4.1), a com- 



parison of the MVFF and Keating- VFF (KVFF) mod- 
els (Sec |4.2[ ), phonon spectrum in Si nanowires (SiNW) 
with free and clamped boundary conditions (Sec. |4.3[ ) 
and lattice thermal conductance using SiNW phonon 
spectrum (Sec 4.4). Conclusions are given in Sec(5} 



2 Theory, Approach and Parameters 

In a given system, the phonons are modeled by solving 
the equations of motion of its atomic vibrations. Since 
VFF is a crystalline model, the dynamical equation for 
each atom 'i' can be written as, 
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where, ARi, Fi and U are the vibration vector of atom 
'i', the total force on atom 'i' in the crystal, and the po- 
tential energy of the crystal, respectively. Equation ([I]) 
indicates that the calculation of the vibrational frequen- 
cies requires a good estimation of the potential energy 
of the system. The next part discusses the calculation 
of U within the MVFF model. 



2.1 Crystal Potential Energy (U) 

The MVFF method JHUHIZ] approximates the potential 
energy U, for a zinc-blende (or diamond) crystal, based 
on the nearby atomic interactions (short-range) [HUH] 
as, 
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where Na, nn(i), and CO Pi represent the total number 
of atoms in one unitcell, the number of nearest neigh- 
bors for atom 'i', and the coplanar atom groups for 
atom 'i', respectively. The first two terms in Eq. ^ are 
from the original KVFF model [7\. The other interac- 
tion terms are needed to accurately describe the phonon 
dispersion in the entire Brillioun Zone (BZ). The terms 
and U^ k represent the elastic energy from bond 
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Fig. 1 The short range interactions used for the calculation 
of phonon dispersion in zinc-blende semiconductors. (a) Bond 
stretching (b) Bond bending (c) cross bond stretching (d) cross 
bond bending-stretching and (e) coplanar bond bending interac- 
tion. 



stretching and bond-bending between the atoms con- 
nected to each other (Fig. 1 a,b). The terms U£l_ bB , 
^bs-bb and U^_ bb represent the cross bond stretch- 
ing [HUH], cross bond bending-stretching [S], and copla- 
nar bond bending [5] interactions, respectively (Fig. [TJ 
c,d,e). The functional dependence of each interaction 
term on the atomic positions is given by, 
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where A9 jik = nj ■ r lk - d ljfi ■ d lkfi , is the angle de- 
viation of the bond between atom 'i' and 'j' and bond 
between atom T and 'k'. The term r^ (d%j o) is the 
non-ideal (ideal) bond vector from atom 'i' to 'j'. The 
coefficients a, ft, 6, 7, and v determine the strength of 
the interactions used in the MVFF model (like spring 
constants). They are used as fitting parameters to re- 
produce the bulk phonon dispersion [8j[9]. The unit of 
these fitting parameters are in force per unit length (like 
Nm^ 1 ). The value of these strength parameters also 
changes according to the deviation of the bond length 
and bond angle from their ideal values. This enables the 
inclusion of the anharmonic properties of the lattice vi- 
brations [50] in this model. Hence, MVFF is sometimes 
referred to as 'quasi-anharmonic' model. 




Fig. 2 (Color online) Three co-planar atom groups (out of 21) 
shown in a bulk zinc-blende unitcell. The groups are (i) 1-2-3-4, 
(ii) 1-2-5-6 and (iii) 1-2-7-8. Atoms 1 and 2 form the bulk unitcell 
used in the calculations. Red (black) atoms are cations (anions). 



Interaction Terms: The primitive bulk unitcell used 
for phonon calculation is made of two atoms (anion- 
cation pair for zinc-blende and 2 similar atoms for di- 
amond). The black dotted box with atom 1 and 2 rep- 
resents the bulk primitive unitcell in Fig. [2] The total 
number of terms in each interaction in Eq. ([3][7| for 
a bulk unitcell are provided in Table [T] Apart from 
the coplanar bond bending interaction (8] all the other 
terms involve nearest neighbor interactions . There are 
21 coplanar (COP) groups present in a bulk zinc-blende 
unitcell which are needed for the calculation of the 
phonon dispersion. For clarity some of these COP groups 
shown using the number combinations in the caption 
of Fig. [2] Each group consists of 4 atom arranged as 
anion(A)-cation(C)-anion(A)-cation(C) (eg. 1(A)-2(C)- 
3(A)-4(C) in Fig. [2]). Details about the coplanar inter- 
action groups are provided in Appendix [A} 



2.2 Dynamical Matrix (DM) 

The dynamical matrix captures the motion of the atoms 
under small restoring force in a given system. In this 
Section we discuss the structure of this matrix. The 
derivation of the DM from the equation of motion is 
given in Appendix |5J The DM calculation is based on 
the harmonic approximation (see Appendix [5]). For the 



Table 1 Number of terms in different interactions of the MVFF 
model in a bulk zinc-blende unitcell (anion-cation pair) 



Interaction Type 


Total terms (anion+cation) 


Bond stretching (bs) 


8 


Bond bending (bb) 


12 


Cross bond stretching (bs-bs) 


12 


Cross bond stretch-bend (bs-bb) 


12 


Coplanar bond bending (bb-bb) 


21 



1 



interaction between two atoms 'i' and 'j', the DM com- 
ponent at atom 'i' is given by, 
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The 9 components of D(ij) are defined as, 
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where Na is the total number of atoms in the unitcell. 
For each atom the size of D(ij) is fixed to 3 x 3. For Na 
atoms in the unitcell the size of the dynamical matrix is 
3Na x 3Na- However, the matrix is mostly sparse. The 
sparsity pattern, fill factor, and other related properties 
of the DM are discussed in Section [3l 

Symmetry considerations in the DM: Under the har- 
monic approximation the dynamical matrix exhibit sym- 
metry properties that can be readily utilized to reduce 
its assembly time. From software development point of 
view this is crucial in optimizing matrix construction 
time, storage and compute times. Due to the continu- 
ous nature of the potential energy U, we have, 
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A closer look at Eq. ( 10 1 shows the following symmetry 
relation, 



D(ij) = D(ji)' Vijtj. 



(11) 



2.3 Boundary conditions (BC) 

To calculate the eigenmodes of the lattice vibration, it 
is important to apply appropriate boundary conditions 
to the DM. In the case of bulk material, the unitcell 
has periodic (Born- Von Karman) boundary conditions 
along all the directions (x,y,z) [SlfTTj since the mate- 
rial is assumed to have infinite extent in each direction. 
However, for nanostructures the boundary conditions 
are different due to the finite extent of the material 
along certain directions. The boundary conditions vary 
depending on the dimensionality of the structure (ID, 
2D or 3D, see Table [5]) for which the dynamical ma- 
trix is constructed. There are 2 types of boundary con- 
ditions; (i) Periodic Boundary condition (PBC) which 
assumes infinite material extent in a particular direc- 
tion and (ii) Finite Edge boundary conditions (FEBC) 
such as open or clamped, which assumes finite material 
extent in a particular direction. Table [2] provides the 
boundary condition details depending on the dimen- 
sionality of the structure used for phonon calculation. 

The use of PBC has been discussed in many papers 
like 8,9,11- In this work we consider the boundary 
conditions associated with geometrically confined nano- 
structures. The vibrations of the surface atoms can vary 
from completely free (free BC) to damped oscillations 
(damped BC). It is shown next that all these cases can 
be handled within one single boundary condition. 

Boundary conditions for nanostructures: The sur- 
face atoms (Fig. [3j hollow atoms) of the nanostructures 
can vibrate in a very different manner compared to 
the inner atoms (Fig. |3j filled atoms) since the surface 
atoms have different number of neighbors and ambient 
environment compared to the inner atoms. The degree 
of freedom of the surface atoms can be represented by 
a direction dependent damping matrix £*, defined in 
Appendix [C] In such a case the dynamical matrix com- 
ponent between atom 'i' and 'j' (D(ij)) is modified to, 



This reduces the total number of calculations re- 
quired to construct the dynamical matrix and speeds up 
the calculations. Also if the matrix is stored for repet- 
itive use, then only one of the symmetry blocks needs 
to be stored. This reduces the memory requirement in 
the software by a factor of 2. Further reduction in the 
construction time of DM can be achieved depending on 
the type of interaction, the symmetry of the crystal, and 
some implementation tricks (not covered in this work, 
see Ref.[2"T] for more discussion). Also the knowledge 
about the underlying symmetry of the matrix can also 
help in the selection of linear algebra approaches which 
can reduce the final solution time (not covered in this 
work) . 



D(ij) = S l D(ij)^ 



(12) 



Table 2 Boundary conditions (BC) in DM based on the dimen- 
sionality of the structure 



Dimensionality 


Periodic BC 


Finite Edge BC 


Bulk (3D) 


3 





Thin Film (2D) 


2 


1 


Wire (ID) 


1 


2 


Quantum Dot (OD) 





3 



5 



Surface atoms 




W(Z) 



Fig. 3 Projected unitcell of a (100) oriented rectangular SiNW 
shown with surface (hollow) and inner (gray filled) atoms. 



2.5 Sound Velocity (V sn d) 

A wealth of information can be extracted from the phonon 
spectrum of solids. One important parameter is the 
group velocity (V grp ) of the acoustic branches of the 
phonon dispersion which gives the velocity of sound 
(V sn d) in the solid. Depending on the acoustic phonon 
branch used for the calculation of V grp , the sound veloc- 
ity can be either, (a) longitudinal (V sn d,i) or (b) trans- 
verse (V sn d,t)- I n solids, V sn d is obtained near the BZ 
center (for q — > 0) where u> ~ q . Thus, V sn d is given 

by, 



V sn d — 



dq 
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where A is the either the transverse or longitudinal po- 
larization of the phonon frequency. 



2.4 Diagonalization of the dynamical matrix 

After setting up the dynamical matrix with appropriate 
BCs the following eigen-value problem must be solved, 



DQ(X,q) =MLu 2 (\,q)Q(\,q), 



(13) 



where, M is the atomic mass matrix. A and q are the 
phonon polarization and momentum vector respectively. 
The term Q(X, q) is a column vector containing all the 
phonon eigen displacement modes u(X, q) associated with 
the polarization A and momentum q. For simplified nu- 



merical calculation slightly modifying Eq. ( 13 ) leads to, 



DQ(X,q)=L0 2 (X,q)Q(X,q) 



(14) 



The detail for obtaining D is outlined in Appendix [d| 
To obtain Eq. (14 1 another step is needed. The time 



dependent vibration of each atom (AR(t)) are repre- 
sented as the linear combination of phonon eigen modes 
of vibration u(X, q) (a complete basis set) as, 



AR l {t) = Y J u P (X, q y^ R ^ 



(15) 



where, P is the size of the basis set and uj the vibration 



frequency of the modes. Using the result of Eq. ( 15 ) in 
the LHS of Eq. yields, 
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After some mathematical manipulations and using Eq. 



( 13 ) we obtain the final eigen value problem given in 



Eq. (14). 



2.6 Thermal Conductance (ct/) 

Another important physical property of the semicon- 
ductors that can be extracted from the phonon spec- 
trum is the lattice thermal conductivity. For a small 
temperature gradient (AT) at the two ends of a semi- 
conductor, the ai is obtained using Landauer approach 
[2"2~] as outlined in Refs. [3"HlfT3"]. For ID nanowires 07 
at a temperature 'T' is given by |4l| 
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where II(w) is the transmission of a phonon branch at 
frequency w, Ti and kg, the reduced Planck's constant 



and Boltzmann constants, respectively. Equation ( 18 ) 



is of general validity and involves a low temperature ap- 
proximation. Scattering causes the conductance to vary 
with the nanowires length (L). In the case of ballistic 
thermal transport the transmission (II(uj)) is always 1 
for all the eigen frequencies. 

2.7 Mode Griineisen Parameter (7.;) 

One of the advantages of using the MVFF model is its 
ability to keep track of the phonon frequency shift under 
crystal stress. Under the action of hydrostatic strain the 
crystal is compressed without changing its symmetry. 
With pressure (P) the phonon frequency shifts, which 
is measured by a unitless parameter called the mode 
'Grzineisen Parameter' given as 
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3 Computational Details 



where, u)^ q is the eigen frequency for the ith branch at q 
momentum point. The terms B, P and V are the volume 
compressibility factor, pressure on the system, and vol- 
ume of the crystal, respectively. Theoretically this pa- 
rameter is extracted by calculating the eigen frequencies 
at ambient condition (P = 0) and at small hydrostatic 
pressure (e = ±0.02) and then taking the difference 
in the calculated frequencies. The modification of the 
force constants under hydrostatic pressure is outlined 
in Ref. [5]. The value of this parameter at high sym- 
metry points (r, X, etc.) in the BZ can be measured 
experimentally by Raman scattering spectroscopy [24]. 

In the remaining Sections we provide the compu- 
tational details, show results on phonon dispersion in 
bulk and nanowires and give some results on V sn d, Ji 
and ballistic <j/ in semiconductor nanowires. 



DM for <100> SiNW, NA = 113 



W = H = 2nm 




VFF-Keating 



Fig. 5 Sparsity pattern of the dynamical matrix used in (a) 
Keating VFF model and (b) MVFF model. SiNW has W = H = 
2nm with 113 atoms in the unitcell. 



This Section provides the computational details involved 
in obtaining the phonon dispersion in semiconductor 
structures. Details about bulk and nanowire (NW) struc- 
tures are provided. 



3.1 Dynamical matrix details 

A primitive bulk zinc-blende unitcell has 2 atoms. This 
fixes the size of the DM for the bulk structure to 6 x 



6 (3Na x 3Na) (see Sec 2.2). However, for the case of 
nanowires, Na varies with shape, size and orientation of 
the wire |13j . In this paper, all the results are for square 
shaped SiNW with (100) orientation. Figure [4] shows 
the variation in Na with the width (W) of Silicon NW 
(SiNW). The number of atoms increase quadratically 
with W. For a 6nm x 6nm SiNW, Na is 1013 which 
means the size of DM is 3,039 x 3,039. Extrapolating 
the Na data gives around 7,128 atoms for a 16nm x 
16nm SiNW resulting in a DM of size 21,384 x 21,384 
(details in Appendix |e| . So, the dynamical matrix size 
increases rapidly with increase in width. 

The increase of the DM size with wire cross-section 
imposes constraint on the structure size which can be 
solved using the atomistic MVFF method. However, the 
entire matrix is quite sparse which can be utilized sig- 
nificantly to expand the physical size of the system that 
can be simulated by using compressed matrix storage 
methods. The qualitative idea about the filling can be 
observed from the sparsity pattern for a 2nm x 2nm 
SiNW dynamical matrix as shown in Fig. [5j The quan- 
titative analysis of the fill fraction of the DM and the 
number of non-zero elements (NZ) in the DM are shown 
in Fig. [6] The non-zero elements in the DM increase 
quadratically with W of SiNW. An estimate for 16nm 
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Table 4 Force constants (N/m) used for phonon dispersion cal- 
culation 



3 4 5 

W = H[nml(log) 



Fig. 7 (a) Time to assemble the DM (time aBrn ) with width (W) 
of (100) oriented SiNW. (b) Time to obtain the eigen solutions 
per k with W for a (100} oriented SiNW DM for 100% and 20% of 
the Eigen spectrum. The timing analysis is done on T6400 Intel 
processor with 2GHz speed. Entire Eigen spectrum along with 
eigen vectors are obtained using the l eig()' function in MATLAB 
1201 . The partial Eigen spectrum is calculated using the l eigs()' 
function in MATLAB [27]. 



x 16nm SiNW gives about 800,117 non-zero elements 
(details in AppendixjE]). However, to get an idea about 
the absolute filling of the DM we define a term called 
the 'fill-factor' given as, 

Total nonzero elements /Size of DM 
NZ 1 

(21) 



fill factor 



(3 x N A ) 

since NZ oc Na (see Appendix Eq. (E.3)) 



Thus, fill factor varies inversely with number of atoms 
in the unitcell. The relation of NZ with NA is provided 
in Appendix [E] 

The percentage fill factor of the DM reduces with 
increasing W of SiNW (Fig. ^. This value is - 0.1% 
for a SiNW with W ~ 25nm (Appendix [E]) . So even 
though the non-zero elements increase with W, DM be- 
comes sparser which allows to store the DM in special 
compressed format like compressed row/column scheme 
(CRS/CCS) [35] enabling better memory utilization. 



3.2 Timing analysis for the computation of DM 

The numerical assembly of the DM takes a consider- 
able time due to the many interactions calculated in the 
MVFF model. The assembly time (time asm ) increases 

Table 3 Resource and timing estimate for larger (100) SiNW 













Time* 


w 


N A 


NZ 


% fill 


tirflCasTn 


per k 


(nm) 






factor 


(sec) 


(hours) 


16 


7128 


800117 


0.423 


238.48 


33.2 


20 


11120 


1252490 


0.224 


370.73 


129.5 


25 


17346 


1.96X10 6 


0.101 


576.91 


505.2 



Material 


Model 


a 


P 


S 


7 


V 


Si 


MVFF [8] 


45.1 


4.89 


1.36 





9.14 


Si 


KVFF [7J 


48.5 


13.8 











Ge 


MVFF [8] 


37.8 


4.24 


0.49 





7.62 



"Time estimates on an Intel T6400, 2GHz processor. 



as Na increases. To give an idea about the timing, the 
dynamical matrix for SiNW with different W are con- 
structed on a single CPU (Intel T6400, 2GHz proces- 
sor). The assembly time is calculated for each width 5 
times to obtain a mean value for the time asm . The er- 
ror bar at each W is the standard deviation from the 
mean time asm (Fig. [7^,). In the present case the assem- 
bly of the DM is done atom by atom which is useful for 
distorted materials as well as alloys. The assembly time 
for the DM in single materials can be reduced dramat- 
ically by the assumption of homogeneous bond lengths 
and a matrix stamping technique [28] . 

After the DM is assembled, it is solved to obtain 
the eigen modes of oscillations. The time needed to 
diagonalize (tdiag) the DM, for each 'k' point, using 
the MATLAB 'eig' function [35] is also calculated (on 
the same processor). The t^ag value varies as the sixth 
power of W as shown in Fig. However, if only 20% 
percent of the Eigen values are calculated the time re- 
quirement now goes by the fifth power (shown by the 
lower line in Fig. [7)3). The Eigen values in this case are 
calculated using the 'eigs' function in matlab [?]. The 
calculation of only 20% of the Eigen spectrum reduces 
the per-k energy calculation time by ~75% for a square 
SiNW with W = H = 6nm. However, the possibility of 
using only the partial spectrum for the evaluation of 
the important lattice parameters (like thermal conduc- 
tance, etc.) is not in the scope of present discussion. For 
the calculation of physical quantities, in this paper, the 
complete Eigen spectrum has been used. 

Extrapolating the data for the computational and 
timing requirement obtained for the smaller SiNWs, 
can provide some estimates about the size and time 
requirement for larger SiNWs (Table [3]). Analytical fits 
for the variation of the size and time parameters with W 
are provided in Appendix [E] The timing requirements 
also help us in the resource requirement estimate for 
a future Bandstructure Lab extension for phonons on 
nanoHUB.org [29] , 



4 Results 

In this Section we show results for the phonon spec- 
trum in bulk and confined semiconductor structures 
using both the MVFF and KVFF models. Also some 
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q [unitless] r ->x->U/K->r ->L 



Fig. 8 Benchmark of simulated bulk phonon dispersion with ex- 
perimental phonon data for (a) Si and (b) Ge. Experimental data 
is obtained using neutron scattering at 80K 30 . 

of the physical properties extracted from the phonon 
dispersions are reported. 

4.1 Experimental Benchmarking 



Table 6 Elastic constants (lO^-Wm -2 ) obtained from the 
MVFF model compared with experimental data |32| for Si and 
Ge. The corresponding errors in the theoretical values are also 
shown. 



Material 


Model 


Cll 


C12 


C44 


Si 


MVFF 


16.80 


6.47 


7.63 


Si 


Expt. 


16.57 


6.39 


7.96 


Error 




~1.4% 


~1.2% 


~4.14% 


Ge 


MVFF 


13.22 


4.84 


6.29 


Ge 


Expt. 


12.40 


4.13 


6.83 


Error 




~6.6% 


~17.2% 


~8% 



point with the experimental data (see Table [7]). The 
MVFF results are in good comparison with the exper- 
imental data. 

An advantage of using a higher order phonon model 
is that both phonon dispersions as well the physical 
parameters can be matched to a good accuracy. Hence, 
MVFF model captures the experimental phonon disper- 
sion as well as the elastic properties in bulk zinc-blende 
material very well. 



The first step to check the correctness of the MVFF 
model is to compare the simulated results with experi- 
mental data. Figure [H] shows the simulated and experi- 
mental [3D] bulk phonon dispersion for (a) Silicon and 
(b) Germanium. The value of the strength parameters 
are provided in Table |4j A very good agreement be- 
tween the experimental and simulated data is obtained. 
To further support the correctness of the MVFF model, 
V sn d is calculated in bulk Si and Ge along (100) direc- 
tion (Table [5]) . The extracted sound velocity compares 
very well with the experimental sound velocity data |31) 
(max error < 10%). 

The comparison of second order elastic constants 
for Si and Ge evaluated using the MVFF model (using 
the formulation provided in Ref. [8] ) with experimental 
data 32 is provided in Table |6j The MVFF derived 
values compare quite well with the experimental data 

E2- 

Another comparison carried out to test the robust- 
ness of MVFF model is the comparison of the mode 
Griineisen parameters at the high symmetry r and X 



Table 5 Sound Velocity in km/sec in Si, Ge bulk and square 
nanowircs with W = H = 6nm. 



Material 


Structure 


V and calc. 


V snd Expt. [31] 


Si 


Bulk Vi [100] 
Bulk V t [100] 
NW Vi 
NW Vt 


9.09 
5.71 
6.51 
4.46 


8.43 (~8%) 
5.84 (~2%) 


Ge 


Bulk Vi [100] 
Bulk V t [100] 
NW V t 
NW Vt 


5.13 
3.36 
3.70 
2.61 


4.87 (~5%) 
3.57 (~6%) 



4.2 Comparison of VFF models 

In this Section we compare the original Keating VFF 
model |7] with the MVFF model to show the need for 
the more elaborate MVFF model. Both computational 
requirement and the physical result comparisons are 
provided in this Section. From computational point of 
view the DM for both the models are quite different 
(Fig. [5]). The difference in the sparsity pattern arises be- 
cause of the coplanar interaction present in the MVFF 
model which takes into account interactions beyond the 
nearest neighbors. The KVFF model has fewer non-zero 
elements compared to the MVFF model. The increase 
in number of NZ elements increases more rapidly in 
the MVFF model vs. the KVFF model (Fig. [9). The 
MVFF model required twice as many matrix elements 
compared to the KVFF model for a 5nm x 5nm SiNW. 
Thus, MVFF model demands more storage space. 

A comparison of the bulk Si phonon dispersion from 
the two models is shown in Fig. 10 The parameters 



for both the models are provided in Table [4] Quali- 
tatively MVFF shows a better match to the experi- 



Table 7 Comparison of the mode Griineisen Parameters for bulk 
Si using the two phonon models. 



7i 


MVFF 


KVFF 


Expt./Abinitio 


Ref. 


1lo,to 


1.05 


0.81 


0.98 ± 0.06 


[2l) 


tTA 


-0.68 


-0.43 


-0.62 


M 


"/La 


0.95 


0.7 


0.85 


[8] 


7lo,la 


1.08 


0.816 


1.03 


1331 




1.25 


0.83 


1.5 ± 0.2 


[241 


Ita 


-1.58 


-0.33 


-1.4 ± 0.3 


[241 



9 



x 10 




Table 8 Comparison of bulk parameters in Si for two models 



Fig. 9 Matrix size and number of non zero elements required by 
the two models. MVFF has more elements needed for accurate 
phonon dispersion. 



mental data compared to the KVFF model. There are 
few important points to be noted in the bulk phonon 
dispersion. The KVFF model reproduces the acoustic 
branches very well near the Brillouin zone (BZ) cen- 
ter but overestimates the values near the zone edge (at 
X and L point in the BZ, Fig. 



10). The MVFF model 



overcomes this shortcoming and reproduces the acous- 
tic branches very well in the entire BZ. Comparison of 
sound velocity along the (100) direction for bulk Si ob- 
tained from both the models show a very good match 
to the experimental data (Table |8| . 

The KVFF model overestimates the optical phonon 
branch frequencies whereas the MVFF model produces 
a very good match to the experimental data (Fig. 10 1 



and Table [8] The comparison of the optical frequency at 
the r point reveals that the KVFF model overshoots 




"0 0.2 0.4 0.6 0.8 f °' 2 °' 4 0-6 0.8 1 

(a) q [unitless]r->X->U/K->r->L (b) 

Fig. 10 Comparison of simulated phonon results with experi- 
mental data (at 80K from 30 ) from the two phonon models (a) 
Keating VFF and (b) Modified VFF. The KVFF model fails to 
reproduce many important features in the experimental data as 
shown by the arrows in (a). The shortcomings are (i) over esti- 
mation of the acoustic mode at X by ~60%, (ii) acoustic branch 
over estimated at L by ~95% and failure to reproduce the correct 
value for the optical branches altogether pointed in (iii) and (iv). 





Model 


\rbulk 
M.lOO 

(km/sec) 


\rbulk 
v t,100 

(km/sec) 


(THz) 


CO 

c= 


MVFF 


9.09 


5.71 


15.49 


CD 


KVFF 


8.35 


5.75 


16.46 


CD 
CD 


Expt. 


8.43 [31] 


5.84 [31] 


15.39 [30] 



the experimental value by ~ 7% whereas the MVFF 
model is higher by only ~ 0.6%. 

The comparison of the mode Griineisen parame- 
ters for bulk Si using the two models is shown in Ta- 
ble [7| KVFF gives wrong values of these parameters 
compared to the experimental values. However, MVFF 
model is able to reproduce the experimental value very 
well. This shows the importance of using a quasi- harmonic 
model to correctly obtain the phonon frequency shifts 
[8,34 . A similar failure of the KVFF model for III-V 
zinc-blende materials have been reported in Rcf. 35 . 

The correct representation of bulk phonons is very 
important since this will affect the phonon spectrum of 
the confined structures. At the same time the physi- 
cal properties like lattice thermal conductivity, phonon 
density of states (DOS), etc. are also affected. Since the 
MVFF model matches the experimental bulk phonon 
data more accurately, though at the expense of addi- 
tional calculations and storage, compared to the orig- 
inal KVFF model, we believe that the MVFF model 
will give better results for phonon dispersion in nano- 
structures. 



4.3 Phonons in nanowires 

After benchmarking the bulk phonon dispersion, the 
same parameters are used to obtain the phonon spec- 



trum in (100) square SiNW (Fig. 111. The result for 
a 2nm x 2nm free-standing SiNW is shown in Fig. 
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a). Some of the key features to notice in the phonon 
dispersion are, (i) presence of two acoustic branches 
(uj(q) ~ q, 1,2 in|TT|a)), (ii) two degenerate modes (3,4 
a)) with oj(q) ~ q 2 , which are called the 'flex- 
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ural modes', typically observed in free-standing nano- 
wires [[6,43,18 and (iii) heavy mixing of the higher en- 
ergy sub-bands leaving no 'proper' optical mode. The 
features are quite different from the bulk phonon spec- 
trum. This will strongly affect other physical properties 
of nanowires extracted using the phonon dispersion. 

In Fig. [Tip , we explore the effect of a substrate on 
which the nanowire may be mounted. Only the bottom 
surface of the SiNW is clamped whereas the other three 
sides have free boundary condition. Using the boundary 
condition method discussed in sec |2.3[ the phonon spec- 
trum in a 2nm x 2nm (100) SiNW are calculated with 
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different damping values (E = 1 (free standing) and 
0.1). The effect of damping is very prominent at the 
BZ edge compared to the zone center. Zone edge fre- 
quencies decrease in energy as the damping increases. 
A reduction of ~ 2.11 x are observed for the zone edge 
frequency for 1st branch at S = 0.1 (Fig. [TT^b)) which 
shows that the NW vibrational energy is decreasing 
more at higher momentum 'q' values. The first four 
branches are very strongly affected, however, the higher 
phonon branches are less affected. 



4.4 Ballistic lattice thermal conductance (&f al ) in 
SiNWs 

The ballistic <7; for square SiNWs is calculated using 
their phonon dispersions. The conductance is calculated 



using Eq. ( 18 ) assuming semi-infinite extensions along 



the wire growth axis (X-axis) and CBC on the periph- 
ery (Y and Z axis) of the wire (Figj3| . Clamping the 
bottom surface affects the a\ al stronger at higher tem- 



perature compared to the lower temperature (Fig. 12 i) 



Figure 12 o shows the a\ al at 300K. The reduction in 



oy from free-standing wire to a clamped wire {E 
0.1 : is ~13.1%. Hence, fixing the surface atoms have 
a strong impact of the lattice thermal conductance in 
SiNW. 



5 Conclusions 

The details for calculating the phonon dispersion in 
zinc-blende semiconductor structures using the modi- 
fied Valence Force Field (MVFF) method have been 




! = 0.1 

8 modes 

2nmX2nm 
(b) SiNW <100> 



0.25 0.5 0.75 1 



q [norm. j 

Fig. 11 (a) Phonon dispersion in (100) oriented SiNW with W = 
H = 2nm. For clarity only the lowest 40 sub-bands are shown, (b) 
Dependence of phonon dispersion on the damping of vibration of 
the bottom surface atoms for S = 1 and 0.1. Reduction of phonon 
energy at the Brillouin zone boundary is stronger compared to 
the zone center. 




100 200 

Temperature [I 



300 



0.4 0.7 1 

Damp Factor (unitless) 



Fig. 12 Ballistic lattice thermal conductivity (er' a! for a 2nm X 
2nm (100) SiNW with different bottom surface damping, af 
drops as damping increases. Inset shows a\ al at 100K, 200K 
and 300K. As the bottom surface changes from free-standing to 
clamped (S = 0.1), a\ al reduces by -12%, ~12.6% and ~13.1% 
at 100K, 200K and 300K, respectively . 



outlined. The MVFF method has been applied to calcu- 
late the phonon spectra in confined nanowire structures 
with varying boundary conditions. The methodology 
and the computational requirements for the method has 
been provided. Comparison of original Keating VFF 
with MVFF shows that, MVFF provides accurate phonon 
dispersion but at the expense of higher computational 
demands. Different VFF models can be used for ob- 
taining the solution for physical quantities depending 
on the type of application, size of the structure and 
the computational resources available. We believe that 
the MVFF model will provide better phonon dispersion 
in ultra-scaled nanostructures. The MVFF method will 
be crucial in understanding and modeling the thermal 
properties of ultra-scaled semiconductor devices. 
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APPENDICES 

A Details of bulk zinc-blende coplanar 
interaction 

The coplanar interactions are important to obtain the flat na- 
ture of the acoustic phonon branches in Si and Ge [8]. There 
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are 21 such interactions in a zinc-blende crystal. The normalized 
locations of all the atoms involved in the coplanar interactions 
are shown in Table [9] The corresponding groups used for bulk 
phonon dispersion calculations are given in Table [To] 



B Derivation of dynamical matrix from the 
equation of motion 

A crystal in equilibrium has zero total force. However, in the pres- 
ence of perturbations like lattice vibrations, etc. a small restoring 
force works on the system. The total force (F tota [) under small 
perturbation is given by the Taylor series expansion as, 

Ftotal = - y, ^ a r, C= Oateqb.) 



1 
2 



dAR, 



d 2 U 



(B.l) 



where, N represents all the atoms present in the system and U 
is the potential energy of the system. In Eq. (JbTTJ first term in 
RHS is zero under equilibrium. The next non-zero term is the 
second term in Eq. |[B7TJ . Under harmonic approximation, only 
the second term is considered and the higher order (anharmonic) 
terms are neglected. Now combining Eq. Q and Eq. ( |B . 1 [ ) one 
can obtain the following, 



3 2 



dt 2 



Ftotal = X/ m 

iejv 

= -- £ 

2 ^— ' 

= DR 



ARi 



d 2 U 
dARidARj 



ARj (B.2) 
(B.3) 

where, D is called the 'Dynamical matrix' and R is a column 
vector of displacement for each atom given as, 



Table 10 Atoms forming the coplanar interaction groups. 4 
atoms in each group. 



No. 


Members 


No. 


Members 


1 


2 


1 


3 


9 


12 


1 


2 


6 


18 


2 


2 


1 


1 


12 


13 


1 


2 


7 


21 


3 


2 


1 


5 


15 


14 


1 


2 


8 


24 


4 


3 


1 


2 


6 


15 


6 


2 


7 


22 


5 


3 


1 


4 


13 


16 


6 


2 


8 


25 


6 


3 


1 


5 


16 


17 


7 


2 


6 


19 


7 


4 


1 


2 


7 


18 


7 


2 


8 


26 


8 


4 


1 


3 


10 


19 


8 


2 


6 


20 


9 


4 


1 


5 


17 


20 


8 


2 


7 


23 


10 


5 


1 


2 


8 


21 


5 


1 


4 


14 


11 


5 


1 


3 


11 













*Atom numbers are same as shown in Table |9[ 



C Treatment of surface atoms 

The damped displacement of the surface atom 'j' can be repre- 
sented the matrix S J given as, 



4 

4 
d 



(C.l) 



Taking into account the individual components the displace- 
ment vector for the atom 'j' we obtain, 



?i = n S [x,y,z] 

This modifies eqn|9]as, 



(C.2) 



(C.3) 



D - 



D(ll) 
D(21) 



D(12) 
D(22) 



D(liV) 
D(2N) 



D(N1) D(N2) ■ ■ ■ D(NN). 



R T = [ARi AR 2 ■ ■ ■ AR N ] 
Definition of D(ij) is given in Eq. (|9j 



(B.4) 



(B.5) 



Combining Eq. | |C.1[ | and | |C.3| the dynamical matrix component 
between atom 'i' and 'j' can be represented as, 



D(ij) 





^x^xy 


f j 
€ il 






A Tjij J 
^y^yx ^x 


i jjij 

y yy 


J 

' y 


t y 1 -'yz 


e? 


i f)ij J 

zx c a? 


i D ij 


t y 


i D 'j 


(? 



(C.4) 



which can be written in a compressed form as, 



Table 9 Normalized atomic coordinates ([x,y, 
used for coplanar interaction calculation. 



[x,y,z 



)/ao) 



No. 


X 


y 


z 


No. 


X 


V 


z 


1* 











14 





0.50 


-0.50 


2* 


0.25 


0.25 


0.25 


15 


-0.50 


-0.50 





3 


0.25 


-0.25 


-0.25 


16 


-0.50 





0.50 


4 


-0.25 


0.25 


-0.25 


17 





-0.50 


0.50 


5 


-0.25 


-0.25 


0.25 


18 


0.25 


0.75 


0.75 


6 





0.50 


0.50 


19 


-0.25 


0.75 


0.25 


7 


0.50 





0.50 


20 


-0.25 


0.25 


0.75 


8 


0.50 


0.50 





21 


0.75 


0.25 


0.75 


9 





-0.50 


-0.50 


22 


0.75 


-0.25 


0.25 


10 


0.50 


-0.50 





23 


0.25 


-0.25 


0.75 


11 


0.50 





-0.50 


24 


0.75 


0.75 


0.25 


12 


-0.50 





-0.50 


25 


0.75 


0.25 


-0.25 


13 


-0.50 


0.50 





26 


0.25 


0.75 


-0.25 



*Belong to the main bulk unitccll used for DM calculation. 



D(ij) = S'Diij)-! 



The value of e x 



(C.5) 

£ [0,1], where completely free surface 



atoms have value 1 and completely tied atoms have value 0. 



D Inclusion of mass in the Dynamical matrix 

In Eq. ( |13| l the mass of the atoms in on the RHS. It is convenient 
to include the mass in DM itself. This modifies the the LHS of 
the equation. The modified DM component between atom 'i' and 
'j' thus, becomes, 



D(ij) 



1 


D a _ 


_L 


1_ 


D a _ 





1_ 


Tjij _ 


_L 


\/mj 


J-^xx 


yrrij 




J-^xy 


Jra j 




L^XZ 


yrrij 


1_ 


D ij - 


1 


1_ 


D a _ 


1 


1 


n« _ 

fyz ^ 


1 




i-Syx 


/™7 


*Jm 


^yy < 






J™j 


1 


Dli - 


_L 


1 


/>.:. 





1 


-DiC 


1 
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Fig. 13 Variation in the number of non-zero (NZ) elements with 
Na for the two phonon models (1) Keating VFF and (2) Modified 
VFF. For both the cases NZ varies linearly with N A . MVFF has 
roughly twice the number of NZ elements compared to KVFF 
model. 



here mi and mj are the masses of atom 'i' and 'j' respectively. 
Eq|D.l|can be written in a compressed manner as, 



D(ij) = M- 1 D(ij)M~ 1 , 
where Mi is given as, 



Mi 



/mi 











m, 

JVTLi 



(D.2) 



(D.3) 



E Fitted analytical expressions for DM 
properties 

— Atoms in a (100) SiNW unitcell: The Na data obtained for 
the square wires till 6nm X 6nm can be fitted to a quadratic 
polynomial given as, 



N A (W) = 27.57W 2 + 4.59W 



(E.l) 



Using Eq. jE.l| l for a 16nm X 16nm SiNW gives around 7128 
atoms. 

Non-zero elements in a (100) SiNW DM: The data for non- 
zero elements in the DM for SiNW with W till 6nm can be 
fitted to a quadratic polynomial given by, 



NZ(W) = 3156W 2 - 495. 5W 



(E.2) 



Using Eq. ]E.2\ for a 16nm X 16nm SiNW yields around 
800117 non-zero elements in the DM. 

Relation of NZ elements to Na in a (100) SiNW: The num- 
ber of non-zero (NZ) elements vary linearly with the number 
of atoms (Fig. |13| |. Fitting the NZ elements with Na for each 
wire under study following relations are obtained, 



NZ KVFF « 109.3 x Na 
NZ MVFF « 205.6 x N A 



(E.3) 
(E.4) 



This shows that the number of NZ elements in MVFF method 
is roughly twice the NZ elements in KVFF model. 



Percentage fill factor for a (100) SiNW DM: The percentage 
fill-factor can be derived using Eq. \E.1\ and \E.2) which 
leads to the following expression, 



%fill - factor(W) i 



0.81 
W* 



oc NA" 



(E.5) 



Eqn |E.5| estimates a 16nm X 16nm SiNW DM is filled only 
0.32%. This shows that the DM matrix is very sparsely filled 
for larger wires. 

Mean assembly time for (100) SiNW DM: The data for mean 
time aS m of the DM for SiNW with W till 6nm can be fitted 
to a quadratic polynomial given by, 



time a sm(W) = 0.9079VK 2 - 0.3789W* 



(E.6) 



Using Eq. ]E.6\ a 16nm X 16nm SiNW DM is estimated to 
be assembled on single CPU in 226.35 sees. 
Eigen solution time per 'k' point for (100) SiNW DM: The 
eigen values are obtained using the 'eig' solver in MATLAB 
|26| . The time needed for the solution of all the eigen values 
(with the eigen vectors) for each momentum vector 'k' can 
be fitted to the following expression, 



time eigen (W) = 5.4 x 10" 3 Vy 6 



(E.7) 



Thus, the solution time goes with the sixth power of W. Also 
this expression gives an estimate time for the solution of one 
k point using the MATLAB eig solver (on a PC) as 1.19e5 
seconds (1.4 days). Thus, for very large systems parallel eigen 
solvers (like SCALAPACK 36 ) as well as finding few eigen 
solutions (using eigs or other sparse eigen solvers, like LA- 
PACK, etc.) is a feasible method. 



References 

1. A. Buin, A. Verma, and M. Anantram, "Carrier-phonon in- 
teraction in small cross-sectional silicon nanowires," Journal 
of Applid Physics, vol. 104, no. 053716, 2008. 

2. A. K. Buin, A. Verma, A. Svizhenko, and M. P. 
Anantram, "Significant Enhancement of Hole Mobility 
in [110] Silicon Nanowires Compared to Electrons and 
Bulk Silicon," Nano Letters, vol. 8, no. 2, pp. 760- 
765, 2008, pMID: 18205425. [Online], Available: |http:| 



|//pubs.acs.org/doi/abs/10.1021/nl0727314| 
N. Mingo and L. Yang, "Phonon transport in nanowires 
coated with an amorphous material: An atomistic Green's 
function approach," Phys. Rev. B, vol. 68, no. 24, p. 245406, 
Dec 2003. 

N. Mingo, L. Yang, D. Li, and A. Majumdar, "Predicting 
the Thermal Conductivity of Si and Ge Nanowires," Nano 
Letters, vol. 3, no. 12, pp. 1713-1716, 2003. 
J. Wang and J.-S. Wang, "Dimensional crossover of thermal 
conductance in nanowires," Applied Physics Letters, vol. 90, 
no. 24, p. 241908, 2007. 

H. Peelaers, B. Partoens, and F. M. Peeters, "Phonon Band 
Structure of Si Nanowires: A Stability Analysis," Nano Let- 
ters, vol. 9, no. 1, pp. 107-111, 2009. 

P. N. Keating, "Effect of Invariance Requirements on the 
Elastic Strain Energy of Crystals with Application to the 
Diamond Structure," Phys. Rev., vol. 145, no. 2, pp. 637- 
645, 1966. 

Z. Sui and I. P. Herman, "Effect of strain on phonons in 
Si, Ge, and Si/Ge heterostructures," Phys. Rev. B, vol. 48, 
no. 24, pp. 17 938-17 953, 1993. 

H. Fu, V. Ozolins, and Z. Alex, "Phonons in GaP quantum 
dots," Phys. Rev. B, vol. 59, no. 4, pp. 2881-2887, 1999. 



13 



10. H. McMurry, A. Solbrig Jr., and J. Boyter, "The use of va- 
lence force potentials in calculating crystal vibrations," Jour- 
nal of Physics and Chemistry of Solids, vol. 28, no. 12, pp. 
2359 - 2368, 1967. 

11. W. Weber, "Adiabatic bond charge model for the phonons 
in diamond, Si, Ge, and ct-Sn," Phys. Rev. B, vol. 15, no. 10, 
pp. 4789-4803, May 1977. 

12. K. Rustagi and W. Weber, "Adiabatic Bond charge model 
for the phonons in A(III)B(V) semiconductors," Solid State 
Communications, vol. 18, pp. 673-675, 1976. 

13. T. Markussen, A. -P. Jauho, and M. Brandbyge, "Heat Con- 
ductance Is Strongly Anisotropic for Pristine Silicon Nano- 
wires," Nano Letters, vol. 8, no. 11, pp. 3771-3775, 2008. 

14. H. L. McMurry, A. W. Solbrig, B. J. K., and C. Noble, "The 
use of valence force potentials in calculating crystal vibra- 
tions," J. Phys. Chem. Solids, vol. 28, pp. 2359-2368, 1967. 

15. J. Zou and A. Balandin, "Phonon heat conduction in a semi- 
conductor nanowire," Journal of Applied Physics, vol. 89, 
no. 5, pp. 2932-2938, 2001. 

16. Y. Zhang, J. X. Cao, Y. Xiao, and X. H. Yan, "Phonon spec- 
trum and specific heat of silicon nanowires," Journal of Ap- 
plied Physics, vol. 102, no. 10, p. 104303, 2007. 

17. X. Li, K. Maute, M. L. Dunn, and R. Yang, "Strain effects 
on the thermal conductivity of nanostructures," Phys. Rev. 
B, vol. 81, no. 24, p. 245318, Jun 2010. 

18. T. Thonhauser and G. D. Mahan, "Phonon modes in Si [111] 
nanowires," Phys. Rev. B, vol. 69, no. 7, p. 075213, Feb 2004. 

19. H. Zhao, Z. Tang, G. Li, and N. R. Aluru, "Quasiharmonic 
models for the calculation of thermodynamic properties of 
crystalline silicon under strain," Journal of Applied Physics, 
vol. 99, no. 6, p. 064314, 2006. 

20. O. L. Lazarcnkova, P. von Allmen, F. Oyafuso, S. Lee, and 
G. Klimeck, "Effect of anharmonicity of the strain energy 
on band offsets in semiconductor nanostructures," Applied 
Physics Letters, vol. 85, no. 18, pp. 4193-4195, 2004. 

21. Z. W. Hendrikse, M. O. Elout, and W. J. A. Maaskant, 
"Computation of the independent elements of the dynamical 
matrix," Computer Physics Communications, vol. 86, no. 3, 
pp. 297 - 311, 1995. 

22. R. Landauer, "Spatial variation of currents and fields due 
to localized scatterers in metallic conduction," IBM J. Res. 
Dev., vol. 1, no. 3, pp. 223-231, 1957. 

23. D. C. Wallace, "Thermodynamics of Crystals," Dover Pub- 
lications, Mineola New York, 1998. 

24. B. A. Weinstein and G. J. Piermarini, "Raman scattering 
and phonon dispersion in Si and GaP at very high pressure," 
Phys. Rev. B, vol. 12, no. 4, pp. 1172-1186, Aug 1975. 

25. J. Dongarra, "Survey of Sparse Matrix Storage Formats," 
1995. [Online], Available: http://www.netlib.org/linalg/ 
html_templates/node90.html 

26. Mathworks, "Matlab erg reference," 2010. [Online]. 
Available: http : / /www. mathworks . com /help/ techdoc / ref / 
|eig.html| 

27. , "Matlab eig reference," 2010. [Online], Availab le: 

http: / /www. mathworks.com/help/tcchdoc/rcf/cigs. html 

28. G. Klimeck, F. Oyafuso, T. B. Boykin, R. C. Bowen, 
and P. von Allmen, "Development of a Nanoelectronic 3-D 
(NEMO 3-D) Simulator for Multimillion Atom Simulations 
and Its Application to Alloyed Quantum Dots," Computer 
Modeling in Engineering and Science ( CMES), vol. 3, no. 5, 
pp. 601-642, 2002. 

29. A. Paul, M. Luisier, N. Neophytou, R. Kim, J. Geng, 
M. McLennan, M. Lundstrom, and G. Klimeck, "Band 
Structure Lab," May 2006. [Online]. Available: |http:| 
/ /nanohub.org/resources/1308 

30. G. Nilsson and G. Nelin, "Study of the Homology between 
Silicon and Germanium by Thermal Neutron Spectrometry," 
Phys. Rev. B, vol. 6, no. 10, pp. 3777-3786, 1972. 



31. "Electronic archive, New Semiconductor Materials - Char- 
acteristics and Properties," Ioffe Physico- Technical Institute 
Website, 2001, http://www.ioffe.ru/SVA/NSM/Semicond/. 

32. O. Madelung, "Semiconductors - HandBook," Springer 3rd 
ed., 2004. 

33. S. de Gironcoli, "Phonons in Si-Ge systems: An ab initio 
interatomic-force-constant approach," Phys. Rev. B, vol. 46, 
no. 4, pp. 2412-2419, Jul 1992. 

34. R. Eryigit and I. P. Herman, "Lattice properties of strained 
GaAs, Si, and Ge using a modified bond-charge model," 
Phys. Rev. B, vol. 53, no. 12, pp. 7775-7784, Mar 1996. 

35. O. L. Lazarenkova, P. von Allmen, F. Oyafuso, S. Lee, and 
G. Klimeck, "An atomistic model for the simulation of acous- 
tic phonons, strain distribution, and Grneisen coefficients in 
zinc-blende semiconductors," Superlattices and Microstruc- 
tures, vol. 34, no. 3-6, pp. 553 - 556, 2003. 

36. L. S. Blackford, J. Choi, A. Cleary, E. D'Azevedo, J. Demmel, 
I. Dhillon, J. Dongarra, S. Hammarling, G. Henry, A. Pe- 
titet, K. Stanley, D. Walker, and R. C. Whaley, ScaLAPACK 
Users' Guide. Philadelphia, PA: Society for Industrial and 
Applied Mathematics, 1997. 



